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Abstract Weak gravitational lensing provides a unique method to map directly the dark matter in the Universe. 
The majority of lensing analyses uses the two-point statistics of the cosmic shear field to constrain the cosmological 
model yielding degeneracies, such as that between erg and £2 m , respectively the r.m.s. of the mass fluctuations 
at a scale of 8Mpc/h and the matter density parameter both at z = 0. However, the two-point statistics only 
measure the Gaussian properties of the field and the weak lensing field is non-Gaussian. It has been shown that the 
estimation of non-Gaussian statistics on weak lensing data can improve the constraints on cosmological parameters. 
In this paper, we systematically compare a wide range of non- Gaussian estimators in order to determine which 
one provides tighter constraints on the cosmological parameters. These statistical methods include skewness, 
kurtosis and the Higher Criticism test in several sparse representations such as wavelet and curvelet; as well as 
the bispectrum, peak counting and a new introduced statistic called Wavelet Peak Counting (WPC). Comparisons 
based on sparse representations show that the wavelet transform is the most sensitive to non-Gaussian cosmological 
structures. It appears also that the best statistic for non-Gaussian characterization in weak lensing mass maps is 
the WPC. Finally, we show that the <T8-f2 m degeneracy could be even better broken if the WPC estimation is 
performed on weak lensing mass maps filtered by the wavelet method, MRLens. 



Weak Lensing, Methods : of its potential to constrain the power spectrum of density 



Key words. Cosmology 
Statistics, Data Analysis 

1. Introduction 



Images distortion measurements of background galax- 
ies caused by large-scale structures provides a direct way 
to study the statistical properties of the growth of struc- 
tures in the Universe. Weak gravitational lensing measures 
the mass and can thus be directly compared to theoreti- 
cal models of structure formation. Most lensing studies use 
the two-point statistics of the cosmic shear field because 
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fluctuations in t he lat e Universe 
2002]; 



Refregier et al 
DahleT 



2005 



e.g 



Maoli et all 



Bacon et all . 120031 : iMassev et al 



2001 



20061 ). Two-point statistics measure the 



Gaussian properties of the field. This is a limited amount 
of information since it is well known that the low redshift 
Universe is highly non-Gaussian on small scales. Indeed, 
gravitational clustering is a non linear process and in 
particular at small scales the mass distribution is highly 
non-Gaussian. Consequently, using only two-point statis- 
tics to set constraints on the cosmological model is lim- 
ited. Better co nstraints can be obtained using 3D weak 



lensing maps (jBernardeau et all Il997t iPen et all 12003 
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Massev et al.1 . 120071 ) . An alternative procedure is to con- 
sider higher-order statistics of the weak lensing shear field 
enabling a characte rization of the non-Gaus s ian nature of 
the signal (see e.g. iTakada and Jainl. 120031: IJarvis et al. , 
2004HKilbinger arid Schneiderl . 120051 : iHamana et alll2004 : 
Donoho and Jinl . l2004 )~ 

In this paper, we systematically compare a range 
of non-Gaussian statistics. For this purpose, we focus 
on the degeneracy between a% and Q m , respectively 
the amplitude of the matter power spectrum and the 
matter density parameter both at z = 0. We investigate 
which statistical method is the most efficient to break 
this degeneracy that exists when only the two-point 
correlation is considered. A wide range of statistical 
methods are systematically applied on a set of simulated 
data in order to characterize the non-Gaussianity present 
in the mass maps due to the growth of structures. Their 
performance to discriminate between different possible 
cosmological models are compared. For the CMB, it has 
been proposed to use statistics such wavelet kurtosis or 
wavelet Higher Criticism to detect clusters and curvelet 
kurtosis or curvelet Higher Cr iticism to d etect anisotropic 
feature such cosmic strings ( Jin et all 20051 ). As weak 
lensing data may contain filamentary structures, we have 
also considered such statistical approaches based on 
sparse representations. 

In section 2, we review the major statistical methods 
used in the literature to constrain cosmological parame- 
ters from weak lensing data. Section 3 describes the simu- 
lations used in this paper, especially how 2D weak lensing 
mass maps of five different models have been derived from 
large statistical samples of 3D N-body simulations of den- 
sity distribution. Section 4 is dedicated to the description 
of our analysis and we introduce the different statistics we 
have studied along with the different multiscale transforms 
investigated. We also present a new statistic that we call 
Wavelet Peak Counting (WPC). In section 5, we present 
our results and finally, section 6 et section 7 present a 
discussion and summaries our conclusions. 

2. Weak lensing statistics and cosmological 
models constraints: state of the art 

2.1. Two-point statistics 

The most common method for constraining cosmol- 
ogy in weak lensing studies uses two-point statistics of 



the shear field calculated either in real or Fourier space. 
In general there is an advantage in using Fourier space 
statistics such as the power spectrum because the modes 
are independent. The power spectrum P K (l) of the 2D con- 
vergence is defined as a function of the modes / by: 



< k{i)k{V) >= {2n) 2 5{i-V)P K (i), 



(1) 



where hat symbols denotes Fourier transforms, <5 is the 
delta function, the brackets denote an average over I. P K {1) 
only depends on / = / and k is the convergence (i.e. which 
is proportional to the projected mass distribution). 

This power spectrum P K (l) can be expressed in terms 
of the 3D matter power spectrum P(k, x) of the mass fluc- 
tuations Sp/ p and of cosmological parameters : 



(fin 



g(x) 



ar(x) 



where a is the expansion parameter, Hq is the Hubble con- 
stant and D. m is the matter density parameter. The corre- 
lation properties are more convenient in Fourier space, but 
for surveys with complicated geometry due to the removal 
of bright stars and artifacts, the missing data need proper 
handling (a review of the di fferent existing methods can be 
found in iPires et all |2008|) . Real space statistics are eas- 
ier to estimate, but require more computational time. The 
following two-point statistics can be related to the under- 
lying 3D matter power spectrum via the 2D convergence 
power spectrum P K (l): 

— The shear variance < 7 2 > : 

An example of real space two-point statistic is the 
shear variance, defined as the variance of the average 
shear 7 evaluated in circular patches of varying radius 
9 S . The shear variance < j 2 > is related to the power 
spectrum P K (l) of the 2D convergence by : 



<r 



>= / > 



(10s) 2 ' 



(3) 



where J n is a Bessel function of order n. This shear 
variance has been used in many weak l ensing analy- 



sis to constrain cosmological parameters (|Maoli et al 
20011 iHoekstra etail 120061 iFu et all l2008fh 



The shear two-point correlation function : 

The shear two-point correlation function is currently 

the most used statistic because it is easy to implement 
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and can be estimated even for complex geometry. It is 
defined as follow : 



=<li{0'h j (0' + 0)>, 



(4) 



where i,j = 1,2 and the averaging is done over pairs of 
galaxies separated by angle = \6\. By parity £ 1-2 = 
£2,1 = and by isotropy £1,1 and £2,2 are functions only 
of 0. The shear two-point correlation functions can be 
related to the 2D convergence power spectrum by : 



£ f (0)=£ M (0)+ 6,2(0) 



dl 
2^ 



IP K (1)J (W), (5) 



These two-point correlation functions are the most 
popular statistical tools and have b een used in the 
most recent weak lensing analysis dBeniamin et a"T 
l2007t iHoekstra etall l2006t iFu et alll200Sh . 



The variance of the aperture mass M ap : 

The aperture mass st atistic has been introduced by 



Schneider et al.l ((1998). It defined a class of statistics 
referred to aperture masses associated with compen- 
sated filters. Several forms of filters have been sug- 
gested which trade locality in real space with local- 
ity in Fou r ier sp ace. Considering the filter defined in 



Schneider (1996) with a cutoff at some scale 6 S . The 



variance of the aperture mass can be expressed as a 
function of the 2D convergence power spectrum by : 



< M 2 ap (6 s ) >= J ^-IP K (1) 



576 Jj(W s ) 
{16.)* ' 



(6) 



This stat istic has been us ed in IVan Waerbeke et al 



2006): IFu et all (120081) . 



(2002); ISem boloni et all (120061 ); IHoekstra et al 



Two-point cosmic shear measurements yield interest- 
ing constraints on the amplitude of the matter power spec- 
trum erg on which the lensing strongly depends. But, when 
deriving joint co nstraints, a degeneracy exists between as 
and fl m (see e.g iMaoli et all l200l l IHoekstra et al.l . l200a 



2.2. Non-Gaussian statistics 

In the standard structure formation model, initial ran- 
dom fluctuations are amplified by non-linear gravitational 
instability to produce a final distribution of mass which is 
highly non-Gaussian. The weak lensing field is thus highly 
non-Gaussian. On small scales, we can observe structures 
like galaxies and clusters of galaxies and on larger scales, 
we observe some filament structures. Detecting these non- 
Gaussian features in weak lensing mass maps can be very 
useful to constrain the cosmological model parameters. 

The three-point correlation function Ci,j,k is the 
lowest-order statistics which can be used to detect non- 
Gaussianity. 



6j,fe(0) =< k(0i)k(0 2 )k(O 3 ) >, 



(7) 



In Fourier space it is called bispectrum and only de- 
pends on distances \l%\ and \ls\: 



BQhlMHsl) ex </s(|li|)ft(IW(l's|)> 



It has been shown 
be obtained with the 



that 



tighter constraints 
point correlation 



(8) 

can 
func- 



tion (Bernardeau et al. 


. 1997; Coorav and Hu. 


Bernardeau et al. . 2003 


Takada and Jainl. I2003L 


Schneider and Lombard] 


, 20031; Schneider et all 


Benabed and Scoccimarrol 20061). 



2001 



2004 



2005 



A simpler quantity than the three-point correla- 
tion function is provided by measuring the third-order 
m oment (skewnes s ) of the smoothed convergence k 



Bernardeau et al 



Jarvis et all 



2004 



19971) or of the aperture mas s M, 



Kilb inger and Schneiderl . [2005l) . 



(If) 



viaon et ai., zuui; noexstra et ai., zuuz; 
20061) . This degeneracy between as and lHamana et all (2 



Another approach to look for non-Gaussianty is to per- 
form a statistical analysis directly on non-Gaussian struc- 
tures like clusters. Galaxy clusters are the largest virial- 
ized cosmological structures in the Universe. They provide 
a unique way to focus on non-Gaussianity present at small 
scales. One interesting statistic is the peak counting that 
searches the number of peaks detected on the 2D con- 
vergence corre s pondi ng to the cluster abundance (see e.g. 



Semboloni et al.l 
ft m is typical of cosmic shear measurements using only 
two-point statistics. 

Two-point statistics are insufficient to characterize 
non-Gaussian features. Non-Gaussianity produced by the 
non-linear evolution of the Universe is of great importance 
for the understanding of the physics of the Universe, and it 
may help to better discriminate the cosmological models. 



2004J). 

The methods to search 



for non-Gaussianity in the 
weak lensing literature mainly concentrate on higher- 
order correlation function. For the CMB, skewness and 
kurtosis of the wavelet coefficients are standard tools as 



Starck et al 



Vielva et al. 2006 



well to check the CMB Gaussiani ty (IVielva et al 
20061; 



20081 ). and it was shown that curvelets (jStarck et al 



2004 



Wiaux et al 
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20031) were useful for the detection o f anisotropic f eature 3.2. 2D Weak lensing mass map 



such cosmic s tring in the CMB (jStarck et all [2004; 



Jin et all 120051 ). 



In this following, we perform a comparison between 
most existing methods in order to find the best higher 
order statistic to constrain cosmological parameters from 
weak lensing data. To explore the effectiveness of a non- 
Gaussian measure we use a battery of N-Body simula- 
tions. By choosing five models whose two-point correla- 
tion statistics are degenerate, we study which statistics 
are able to distinguish between these models. 

3. Simulations of weak lensing mass maps 

3.1. 3D N-body cosmological simulations 

We have run realistic simulated convergence mass 
maps derived from N- body cosmolog ical simulations using 



20021 ). The cosmological 



the RAMSES code (|Tevssieii 
models are taken to be in concordance with the ACDM 
model. We have limited the model parameters to a 
realistic range (see Table [IJ choosing five models along 
the (cr 8 , fi m )-degeneracy discussed in ^2.11 



Model 


Box 




n L 


h 


era 


model 1 


165.8 


0.23 


0.77 


0.594 


1 


model 2 


159.5 


0.3 


0.7 


0.7 


0.9 


model 3 


152.8 


0.36 


0.64 


0.798 


0.8 


model 4 


145.7 


0.47 


0.53 


0.894 


0.7 


model 5 


137.5 


0.64 


0.36 


0.982 


0.6 



Table 1. Parameters of the five cosmological models that 
have been chosen along the (ag, fi m )-degeneracy. The box 
side is in Mpc/h, the simulations have 256 3 particles, Vt m 
is the matter density parameter, is the dark energy 
density parameter, h is equal to iJo/100 where H$ is the 
Hubble constant and ag the amplitude of the matter power 
spectrum. 



In N-body simulations, which are commonly used in 
cosmology, the dark matter distribution is represented us- 
ing discrete massive particles. The simplest way to treat 
these particles is to map their positions onto a pixelised 
grid. In the case of multiple sheet weak lensing, we do this 
by taking slices through the 3D simulations. These slices 
are then projected into 2D mass sheets. 

The effective convergence can subsequently be calcu- 
lated by stacking a set of these 2D mass sheets along the 
line of sight, using the lensing efficiency function. This is 
a proc edure that has been used before by I Vale and White 
(|2003h . where the effective 2D mass distribution n e is cal- 
culated by integrating the density fluctuation along the 
line of sight. Using the Born approximation which neglects 
the facts that the light rays don't follow straight lines, the 
convergence can be numerically expressed by : 



2^2 



Xi) 



Xoa(Xi 



n p R 2 



A 



r fi 



(9) 



where Hq is the Hubble constant, fl m is the density of mat- 
ter, c is the speed of light, L is the length of the box, x 
are co-moving distances, with xo being the co-moving dis- 
tance to the source galaxies. The summation is performed 
over the i th box. The number of particles associated with 
a pixel of the simulation is n p , the total number of par- 
ticles within a simulation is N t and s = L p / L where L p 
is the length of the plane doing the lensing. R is the size 
of the 2D maps and Ar/$ = r2 ^ ri where r\ and are 
co-moving distances. 

Using the previous 3D N-body simulations, we have 
derived 100 realizations for the five models. Fig. Q] shows 
one realization of the convergence for each of the 5 models. 
In all cases, the field is 3.95° x 3.95°, downsampled to 
1024 2 pixels and we assume that the sources lie at exactly 
z = 1. At large scale, the map shows clearly a Gaussian 
signal. On the contrary, at small scales, the signal is clearly 
dominated by clumpy structures (dark matter halos) and 
is therefore highly non- Gaussian. 



For each of our five models, we run 21 N-Body simula- 
tions, each containing 256 3 particles. We refined the base 
grid of 256 3 cells when the local particle number exceeds 
10. We further refined similarly each additional levels up 
to a maximum level of refinement of 6, corresponding to 
a spatial resolution of 10 kpc.h -1 . 



3.3. 2D Weak lensing noisy mass map 

In practice, the observed shear 7, is obtained by av- 
eraging over a finite number of galaxies and is therefore 
noisy. The noise arises both from the measurement errors 
and from the intrinsic ellipticities dispersion of galaxies. 



0.2 



0.4 0.6 0.8 1 



■ * 

•..(.•v. 




Figure 1. Upper left, the 5 cosmological models along the (as, ri m )-degeneracy. Upper right, one realization of the 
convergence k for model 1 (erg = 1 and Q m — 0.23); middle left, for model 2 (as = 0.9 and fl m = 0.3); middle right, 
for model 3 (erg = 0.8 and fi m = 0.36); bottom left, for model 4 (<jg = 0.7 and fl m = 0.47) and bottom right for model 
5 (as = 0.6 and Q m = 0.64). Each map is 1.975 x 1.975 degrees down-sampled to 512 x 512 pixels 
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As a good approximation, we modeled the noise as an un- 
corrected Gaussian random field with variance : 



(10) 



.4. 



where A is the pixel size in arcmin 2 , n g the average num- 
ber of galaxies per arcmin 2 and a e the rms of the shear 
distribution. We assume a e ~ 0.3 per shear component 
and n g — 100 gal/arcmin 2 as is ap proximately found fo r 



space-based weak lensing surveys (jMassey et all 120041 ) . 



An estimate of the noisy convergence K n is obtained usin g 
the least square estimator defined in Starck et al. (2006). 



4. Cosmological model discrimination framework 

In this study, in order to find the best statistic 
to break the (as, f2 m )-degeneracy, we are interested 
in comparing different statistics estimated in dif- 
ferent representations using the set of simulated 
data described in the previous section. 

For this purpose, for each statistic, we need to 
characterize, the discrimination obtained between 
each couple of models. The best statistic will be 
the one that maximizes the discrimination for all 
the couples of models. 

4.1. Characterization of the discrimination 

To find the best statistic, we need to charac- 
terize quantitatively for each statistic the discrim- 
ination between two different models mj and m-x. 
One way to proceed is to define a discrimination 
efficiency that expresses the ability of a statistic 
to discriminate in percentage. Then, we need to 
define for each individual statistic, two different 
thresholds (see Fig. [2]): 

- a threshold T\ : if the estimation of a given statis- 
tic in a map k% is below T\ the map n% is classified 
belonging to the model mi and not if it is above 

- a threshold T2 : if the estimation of a given 
statistic in a map K\ is above Ti the map /t, 
is classified belonging to the model m.2 and not 
if it is below. We have used a statistical tool 
called FDR (False Discov e ry Ra te) introduced by 
Beniam ini and Hochberd ( 1995h to set these two 
thresholds (T\ and T2) correctly (see Appendix A). 




Figure 2. The following two distributions corre- 
spond to the histogram of the values of a given 
statistic estimated on the 100 realizations of model 
1 (mi) and on the 100 realizations of model 2 (77x2). 
The discrimination achieved with this statistic be- 
tween mi and m,2 is rather good : the two distribu- 
tions barely overlap. To characterize more quanti- 
tatively the discrimination, the FDR method has 
been used to estimate the thresholds T\ and T2. 
A false discovery rate (a) equal to 0.05 has been 
chosen. Then a discrimination efficiency can be de- 
rived. 



This FDR method is a competitive tool to set 
a threshold in an adaptive way without any as- 
sumption, given a false discovery rate (a). The 
false discovery rate is given by the proportion of 
false detections over the total number of detec- 
tions. The threshold being estimated, we can de- 
rive a discrimination efficiency for each statistic. 
The discrimination efficiency measures the abil- 
ity of a statistic to discriminate a model from an- 
other one by calculating the percentage of detec- 
tions (true or false) over the total number of sam- 
ples. It corresponds basically to the part of the 
distribution that doesn't overlap. The more the 
distributions overlap, the lower the discrimination 
efficiency will be. 

The Fig. [2] represents the dispersion of the val- 
ues of a given statistic estimated on the 100 real- 
izations of the model 1 (on the left) and on the 100 
realizations of the model 2 (on the right). The two 
distributions barely overlap, it indicates a good 
discrimination that is to say the two models can 
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easily be separated with this statistic. To be more 
quantitative, a threshold has been set for each dis- 
tribution to estimate a discrimination efficiency 
that corresponds to the part of the distribution 
delimited by the hatched area. 

The formalism of the FDR method ensures that 
the yellow area delimited by T\ (resp. T%) that cor- 
responds to the false detections will be small. 

4.2. A set of statistical tools 

The first objective of this study is to compare different 
statistics to find the one that sets tighter constraints on 
the cosmological parameters. The two-point statistics that 
contain all the information about Gaussian signal lead- 
ing to the (as, ri m )-degeneracy, we have opted for statis- 
tics currently used to detect non-Gaussianity in order to 
probe the non-linear process of gravitational clustering. 
The statistics that we have selected are the following : 

— Skewness (S K ): 

The skewness is the third-order moment of the conver- 
gence k and is a measure of the asymmetry of a distri- 
bution.The convergence skewness is primarily due to 
rare and massive dark matter halos. The distribution 
will be more or less skewed positively depending on the 
abundance of rare and massive halos. 

— Kurtosis (K K ): 

The kurtosis is the fourth-order moment of the conver- 
gence k and is a measure of the peakedness of a dis- 
tribution. A high kurtosis distribution has a sharper 
"peak" and flatter "tails", while a low kurtosis distri- 
bution has a more rounded peak with wider "shoul- 
ders" . 

— Bispectrum (B K ): 

Recently, there has been a lot of theoretical works on 
the three-point correlation function to constrain the 
cosmological parameters. But the direct computation 
of three-point correlation function takes too long for 
our large maps. We have then used its Fourier analog: 
the Bispectrum that has been introduced §2.21 And we 
have consider the equilateral configuration. 

— Higher Criticism (HC): 

HC is a recent ly deve loped statistic, proposed by 
Donoho and Jinl ( 2004 ). It is a measure of non- 



large HC value implies non-Gaussianity. We consider 
the two different forms of HC (see Appendix B) : 

- HC* 

- HC+ 

— Peak counting (P c ): 

We now investigate the possibility of using the peak 
counting to differentiate among cosmological models. 
By peak counting (or cluster count), we mean the num- 
ber of halos that we can detect per unit area of the 
sky (identified as peak above a mass threshold in mass 
maps) . This cluster count enables to constrain the mat- 
ter pow er spectrum normalizat ion as for a given f2 m 
(see e.g. lBahcall and Fanlll998l ) and the formalism ex- 
ist to predict t he peak count to a g iven cosmological 
model (see e.g. Hamana et all [2004 ) . 

- Wavelet Peak Counting (WPC): 

We introduce a new statistic that we call Wavelet Peak 
Counting (WPC). This statistic consists in estimating 
a cluster count per scale of a wavelet transform. It 
means we have roughly a cluster count depending on 
the size of the clusters. We show in the following that 
the WPC is better than the peak counting to charac- 
terize the non-linear structure formation process. 

All the statistics are estimated on the 2D convergence 
k map. In practice the complex geometry of surveys gives 
the convergence n inferred from the shear field ji the same 
complex geometry. But it beco mes possibl e with the in- 
painting method developed in IPires et al.1 (|2008h to re- 
construct a full convergence k map. The solution that is 
proposed enables to fill-in judiciously masked regions so 
as to reduce the impact of missing data on the estima- 
tion of the power spectrum and of higher-order statistical 
measures while requiring 0(N log N) operations. 

4.3. Representations 

The second objective of the present study is to com- 
pare different transforms to find the sparsest representa- 
tion of weak lensing data that makes the discrimination 
easier. Some studies on CMB data have used multiscale 



Gaussianity obtained by looking at the maximum de- 
viation comparing the sorted p-values of a distribution 
with the sorted p-values of a normal distribution. A 



methods to detect non-G aussianity (jAghanim and Fornil . 
19991: IStarck et all 120041) . Weak lensing maps exhibit both 
isotropic and anisotropic features. These kind of features 
can be better represented in some basis functions. Indeed, 
a transform is optimal to detect structures which have 
the same shape of its basis elements. We have thus tested 
different representations : 
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— The Fourier transform 

— The anisotropic bi-orthogonal wavelet transform 

We expect the bi-orthogonal wavelet transform to be 
optimal for detecting mildly anisotropic features. 

— The isotropic " a trous" wavelet transform 

This wavelet transform is well adapted to the detec- 
tion of isotropic features such as the clumpy structures 
(clusters) of the weak lensing data. 

— The ridgelet transform 

The ridgelet transform has been developed to process 
images including ridge elements, and so provides a 
good representation of perfectly straight edges. 

— The curvelet transform 

The curvelet transform allows to approximate curved 
singularities with few coefficients and then provides a 
good representation of curved structures. 

A description of these transforms is available in Appendix 
C. We have to notice that the estimation of a given statis- 
tic at each resolution level of a representation can be con- 
sidered as an individual statistic by itself. 

5. Analysis and results 

5.1. Treatment of the noise 

As explained previously, the weak lensing mass maps 
are measured from the distortions of a finite number of 
background galaxies and therefore suffer from shot noise. 
Furthermore each galaxy provides only a noisy estimator 
of the distortion field. We have added the expected level 
of Gaussian noise (see ^3.3|) to simulations of weak lensing 
mass maps to obtain simulated noisy mass maps corre- 
sponding to space observations. Fig. [3] shows a noiseless 
simulated mass map (left) and a noisy simulated mass map 
(right) corresponding to space observations. 

The noise has an impact on the estimated statistics 
and therefore needs to be considered. 



5.1.1. No filtering 

We begin by applying the different transformations di- 
rectly to the noisy data and for each representation, we 
have estimated the statistics described in the previous sec- 
tion, with the exception of cluster count and WPC which 
required filtering. 

We expect that the noise will make the third and 
fourth order statistics tend to zero. Indeed the more noisy 



the data are, the more the distribution will look like a 
Gausssian. 

5.1.2. Gaussian filtering 

As a second step, we have performed a Gaussian filter- 
ing that is obtained by convolving the noisy mass maps 
K n with a Gaussian window G with a standard deviation 



Kq = G * K r , 



(11) 



Fig.|4]shows on the left, the original map without noise 
and on the right, the result obtained by Gaussian filtering 
of the noisy mass map displayed in Fig. [3] (right). The 
quality of the filtering depends strongly on the value of 
ere?. For the simulations, an optimal value is about 0.9 
arcmin). We then computed all our statistics on this new 
set of filtered mass maps. 

5.1.3. MRLens filtering 

Finally, we used a method of non-linear filtering based 
on the wa yelet representation : the MRLens filtering pro- 
posed bv IStarck et al.l (|2006l) . The MRLens filtering 
is based on the Bayesian ideas that provides the 
means to incorporate prior knowledge in data anal- 
ysis. Choosing the prior is one of the most criti- 
cal aspects of Bayesian analysis. Here a multiscale 
entropy prior is used. A full description of the 
method is given in the Appendix D. The MRLens 
software that we used is available at the following 
address : http: / / ww w-irm.cea.fr/Ast/8 78.html. 

It was shown in Starck et al. (2006) that this 
method outperforms several standard techniques 
(Gaussian filtering, Wiener filtering, MEM filter- 
ing, etc.). Fig. [5] shows on the left, the original map 
without noise and on the right, the result of a 
FDR filtering on the noisy mass map displayed Fig. 
[3] (right). The visual aspect indicates that many 
clusters are reconstructed and the intensity of the 
peaks are well recovered. 

As done before, all the statistics have been es- 
timated on these MRLens filtered mass maps. By 
reconstructing essentially the clusters, we can anticipate 
that this MRLens method will help statistics like peak 
counting or WPC more than statistics that focus on the 
background. 





Figure 3. Left: noiseless simulated mass map and right: simulated noisy mass map that we should obtain in space 
observations. The field is 1.975° x 1.975°. 




Figure 4. Left, noiseless simulated mass map, and right, filtered mass map by convolution with a Gaussian kernel. 
The field is 1.975° x 1.975°. 



5.2. Results els out of the full set of 5 models. This values are 

given Table [H Table [5] and Table [7] for three dif- 
5.2.1. The discrimination methodology ferent statistics. 

A mean discrimination efficiency for each in- 
As explained section 14.11 for each statistic de- dividual statistic can be estimated by averaging 
scribed in the previous section, we can derive a the discrimination efficiency across all the pairs of 
discrimination efficiency between each two mod- models. For statistics estimated in multiscale rep- 
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Figure 5. Left, noiseless simulated mass map, and right, filtered mass map by the FDR multiscale entropy filtering 
The field is 1.975° x 1.975°. 



resentations, the mean discrimination efficiency is 
calculated for each scale and it is the better one 
that is considered. 

Tables |2l [4] and [6] show the mean discrimina- 
tion efficiency reached by a given statistic and 
a given transform, respectively for (i) unfiltered 
mass maps, (ii) Gaussian filtered mass maps and 
(iii) MRLens mass maps. 

The mean discrimination efficiency of the Table 
[3] is about 40% and corresponds to the value at po- 
sition (1,3) of Table [2] (i.e. the skewness of wavelet 
coefficients). And, the mean discrimination effi- 
ciency of Table [7] is about 97% and corresponds 
to the value at position (3,5) of Table [6J 

A full discrimination between the five models 
is obtained if the mean discrimination efficiency is 
about 100%. If it is between 80% and 100% the dis- 
crimination is possible except for adjacent models. 
Less than 40% there is no discrimination possible 
even for distant models. 

5.2.2. Discrimination in noisy mass maps 

The mean discrimination efficiency obtained for unfil- 
tered mass maps is displayed in Table [21 The peak count- 
ing and WPC differ from the others. They can not be 
computed on unfiltered mass maps because the clusters 
can not be extracted from noisy mass maps. Another point 



is that the bispectrum by definition can only be estimated 
in the Fourier domain. 

Without filtering the results are poor and no discrimi- 
nation can be achieved in Direct space. Indeed, the Signal- 
to-Noise Ratio is weak as can be seen in Fig. [3] (right) . The 
non-Gaussian signal is hidden by the Gaussian noise. 

The different transforms appear to be inefficient 
at bring out the non-Gaussianity features, except, the 
Isotropic Wavelet Transform (see Table [2]) which performs 
rather well, whatever the statistic. This is likely because 
it is optimal to detect clusters. Indeed, the clusters are a 
direct probe of non-Gaussianity and by concentrating the 
cluster information, the isotropic wavelet transform seems 
to be a better representation for non-Gaussianity. 

The skewness in the Wavelet transform representa- 
tion appears to be the best statistic in unfiltered mass 
maps. Table [3] shows the discrimination efficiency ob- 
tained with the skewness at the second scale of an 
Isotropic Wavelet Transform being the scale that reaches 
the best discrimination. We can see that the discrimina- 
tion is only achieved between the farthest models, which 
is quite poor. It is illustrated on Fig. [6] where you 
can see that the 5 distributions are overlapping. 
Some groups have already used the skewness aperture 
mas s statistic to try to break the cr« - n m degeneracy, (se e 



e.g. iKilbinger and Schneiderl . l200a Ijarvis et all 12004! ) . 



This processing consists of convolving the noisy signal by 
Gaussian windows with different scales and is quite simi- 
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Basis/Statistics 


S 


K 




HCZ 




Pc 


WPC 


Direct space 


21.0 


1.6 


1.8 


1.9 


X 


X 


X 


Fourier space 


0.7 


10.7 


3.0 


21.3 


5.75 


X 


X 


Isotropic Wavelet Transform 


40.5 


29.3 


24.9 


34.0 


X 


X 


X 


Bi-orthogonal Wavelet Transform 


7.3 


10.8 


4.1 


2.9 


X 


X 


X 


Ridgelet Transform (block = 16) 


6.2 


3.7 


7.2 


5.6 


X 


X 


X 


Ridgelet Transform (block = 32) 


8.3 


15.2 


9.5 


20.4 


X 


X 


X 


Fast Curvelet Transform 


1.1 


4.2 


2.6 


4.2 


X 


X 


X 



Table 2. Mean discrimination efficiencies (in percent) achieved on noisy mass maps with a False Discovery Rate 
a = 0.05. 



lar to an Isotropic Wavelet Transform. They showed that 
by combining the second and third-order statistics, the 
degeneracy can be diminished but not broken. 



m odel 
rn odel 
rri cxJc;! 
r~r~i odel 
m odel 



IT. 




TL 



HI 



rn~i 



ru 



n^tn;-r nW. n n 



Figure 6. Distribution of the skewness calculated 
from the second scale of an Isotropic Wavelet 
Transform on the simulated realizations of the 5 
models. It illustrates the results of the Table 3. 
No discrimination is possible except between the 
farthest models (i.e. between model 1 and model 
5) 



ilar. Some statistics are a bit improved by the noise re- 
moval while others become worse. 

By contrast, the peak counting and WPC that can 
now be estimated on these filtered mass maps perform 
well. Table [5] shows the discrimination efficiency obtained 
with peak counting estimated on Gaussian filtered mass 
maps. We can see that except for adjacent models, dis- 
crimination is now possible. We can verify this results 
by looking at the 5 distributions displayed Fig. \7\ 
Indeed, the distributions barely overlap for no ad- 
jacent models. 

The ability of weak lensing cluster count (peak count- 
ing) to discriminate between the five different models cho- 
sen along the degeneracy can be explained. Considering 
that the dark matter lies at z ~ 0.5, which is the maxi- 
mum lensing efficiency for background galaxies at z = 1, 
and assuming a constant dark energy model, the number 
density of massive clusters in the weak lensing mass maps 
is sensitive to both the amplitude of the mass fluctuations 
as and the matter density parameter Q m , both at z = 0. 
If instead of erg, we consider erf- 5 that corresponds to 
the amplitude of the fluctuations of the projected weak 
lensing mass map. The cr|- 5 is now a constant along 
the five models because the five corresponding weak lens- 
ing power spectrum are undistinguishable. This leaves the 



5.2.3. Discrimination in Gaussian filtered mass maps pm parameter that drives the structure formation (see e.g 



To increase the Signal-to-Noise Ratio, we have applied 
a Gaussian filtering to the noisy simulated mass maps. 
Table H] shows the results. 

After Gaussian filtering, the noise is removed but the 
structures are over-smoothed. Except for the direct space 
where the results are clearly improved by the noise re- 
moval, the results after a Gaussian filtering are quite sim- 



Bahcall and FarJ . 1998h . A small fl m makes the structures 
form earlier and a large Cl m makes the structures form 
later. Then, with a small f2 m , the abundance of massive 
clusters at z ~ 0.5 are more significant (see Fig. Q] up- 
per right) than for a large f2 m The cluster count can then 
be used to discriminate cosmological model. The massive 
clus ter abundance ha s alrea dy been used to probe fl m (see 
e.g. iBahcall and Fanl . [l998h . 
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model 1 


model 2 


model 3 


model 4 


model 5 


model 1 


X 


4 


55 


86 


99 


model 2 


6 


X 


8 


52 


88 


model 3 


23 


3 


X 


6 


77 


model 4 


81 


7 


2 


X 


13 


model 5 


99 


81 


12 


9 


X 



Table 3. Discrimination efficiencies (in percent) achieved on uiffiltered mass maps with the skewness estimated at the 
second scale of an Isotropic Wavelet Transform. 



Basis/Statistics 


S 


K 






Bi 


Pc 


WPC 


Direct space 


39.2 


34.3 


28.8 


40.2 


X 


79.3 


X 


Fourier space 


2.2 


2.2 


3.8 


2.9 


6.75 


X 


X 


Isotropic Wavelet Transform 


31.9 


29.5 


24. 


38.4 


X 


X 


79.2 


Bi-orthogonal Wavelet Transform 


12.2 


10.5 


5.1 


6.8 


X 


X 


X 


Ridgelet Transform (block = 16) 


7.4 


20.2 


22.7 


35.3 


X 


X 


X 


Ridgelet Transform (block = 32) 


2.5 


7.3 


4. 


9.1 


X 


X 


X 


Fast Curvelet Transform 


1.3 


4.7 


4.8 


8.3 


X 


X 


X 



Table 4. Mean discrimination efficiencies (in percent) achieved on Gaussian filtered mass maps with a False Discovery 
Rate a = 0.05. 





model 1 


model 2 


model 3 


model 4 


model 5 


model 1 


X 


50 


95 


100 


100 


model 2 


46 


X 


43 


96 


100 


model 3 


94 


8. 


X 


63. 


100 


model 4 


100 


97 


50 


X 


78 


model 5 


100 


100 


99 


63 


X 



Table 5. Discrimination efficiencies (in percent) achieved on Gaussian filtered mass maps with the peak counting 
statistic on direct space given a False Discovery Rate a = 0.05. 



5.2.4. Discrimination in MRLens filtered mass maps 

Table [5] shows the results obtained in the case where 
the MRLens filtering scheme is applied to the noisy simu- 
lated mass maps. 

After MRLens filtering, the sensitivity of all trans- 
forms is greatly improved. But, the Isotropic Wavelet 
Transform remains the best transform certainly helped by 
the MRLens filtering that uses this transform. 

The best statistics remains the peak counting, which 
is also helped by the MRLens filtering that reconstructs 
essentially the clusters. But the others statistics also 
achieve quite good results in these MRLens filtered mass 
maps compared with the one with Gaussian filtered mass 



maps. This is likely because the MRLens filtering by 
favoring the clusters reconstruction helps all the statistics 
that look for non-Gaussianity. 

The best result is obtained with WPC on the third 
scale of the Isotropic wavelet transform see Table[7] Fig. [8] 
shows the 5 distributions that barely overlap. This 
statistic allows us to discriminate between different models 
even for adjacent models for which the discrimination is 
challenging. 

The comparison of these results with the results ob- 
tained on noisy mass maps by the skewness in a wavelet 
representation (Tabled show that the accuracy on ct 8 and 
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Basis/Statistics 


S 


K 




HCZ 


Bi 


Pc 


WPC 


Direct space 


65.2 


54.5 


72.9 


8.4 


X 


93.2 


X 


Fourier space 


3.2 


3.4 


30.0 


38.4 


9.95 


X 


X 


Isotropic Wavelet Transform 


81.7 


68.1 


74.8 


74.1 


X 


X 


97.3 


Bi-orthogonal Wavelet Transform 


52.2 


67.3 


29.7 


65.0 


X 


X 


X 


Ridgelet Transform (block = 16) 


45.1 


54.5 


77.1 


37.2 


X 


X 


X 


Ridgelet Transform (block = 32) 


61.5 


52.6 


70.3 


69.0 


X 


X 


X 


Fast Curvelet Transform 


25. 


65.5 


68.8 


63.3 


X 


X 


X 



Table 6. Mean discrimination efficiencies (in percent) achieved on MRLens filtered mass maps with a False Discovery 
Rate a = 0.05. 





model 1 


model 2 


model 3 


model 4 


model 5 


model 1 


X 


86 


100 


100 


100 


model 2 


87 


X 


94 


100 


100 


model 3 


100 


92 


X 


94 


100 


model 4 


100 


100 


93 


X 


99 


model 5 


100 


100 


100 


100 


X 



Table 7. Discrimination efficiencies (in percent) achieved on MRLens filtered mass maps with WPC at the third scale 
of an Isotropic Wavelet Transform. 
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Figure 7. Distribution of the peak counting esti- 
mated directly on the simulated realizations of the 
5 models. It illustrates the results of the Table 5. 
The discrimination is possible except between ad- 
jacent models (that is to say between model 1 and 
model 2, model 2 and model 3, model 3 and model 
4, model 4 and model 5). 



J7 m is greatly improved using WPC estimated on MRLens 
filtered mass maps. 



Figure 8. Distribution of the Wavelet Peak 
Counting estimated at the third scale of an 
Isotropic Wavelet Transform on the simulated re- 
alizations of the 5 models. It illustrates the results 
of the Table 7. We obtain a good discrimination 
even for adjacent models. 

6. Discussion 

As stated earlier, the formalism of the halo model 
provides a prediction for the number of clusters con- 
tained in a given fie l d for a given cosmologica l model 
(jPress and Schechterl . 11974 ISheth and Tormenl . Il999t 
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Hamana et all 120041 ). However, we have to consider that 
just a fraction of the clusters, present in the sky, will be 
detected. It follows that we have to take into account the 
selection effects coming from the observation quality and 
the data processing method. The solution that is currently 
used consists in modeling the selection effects by estimat- 
ing the selection function. An analytic model can be done 
by considering all the selection effects. An alternative con- 
sists in using a Monte Carlo approach that enables us to 
take into account the entirety of the selection effects that 
could not be considered in the analytic approach. This 
study will be done in a future work. 

The selection function being specified, the connection 
with observations and theory is straightforward. The cos- 
mological parameters can thus be estimated from WPC ' . 

Having a perfect discrimination between the 5 cosmo- 
logical models with WPC upper limits error on the cosmo- 
logical parameters can be given by considering the spacing 
between two adjacent models. For space observations cov- 
ering 4 square degrees, the upper limit error in as is 8%, 
in the range of as € [0.6, 1]. And the upper limit error in 
ft m is 12%, in the range of Sl m £ [0.23,0.64]. In future 
work, an accurate estimation of the error should be done. 

7. Conclusion 

Using only two-point statistics to constrain the cosmo- 
logical model yield to various degeneracies between cos- 
mological parameters such as the (o"s-f2 m )-degeneracy. In 
this paper, we survey a range of non-Gaussian statistics 
to set tighter constraints on cosmological parameters. For 
this purpose, we have run N-body simulations of 5 mod- 
els selected along the (o-g-Q m )-degeneracy and we have 
examined different non-Gaussian statistical tools in differ- 
ent representations in order to compare the discrimination 
power of each one. Using non-Gaussian statistics, we have 
looked for non-Gaussian signal present at small scales that 
is due to gravitational clustering. 

The main conclusions of our analysis are the following: 

1. The isotropic wavelet transform has been found to be 
the best representation of the non-Gaussian structures 
in weak lensing data allowing a better discrimination. 

2. We have shown that a wavelet transform denoising 
method like the MRLens filtering, which reconstructs 
essentially the non-Gaussian structures (the clusters), 
helps the statistics to better characterize the non- 
Gaussianity. 



3. We have introduced a new statistic called Wavelet 
Peak Counting (WPC) that consists in estimating a 
cluster count per scale of an isotropic wavelet trans- 
form. 

4. WPC has been found to be the best statistic compared 
to the others that we have tested (skewness, kurto- 
sis, bispectrum, HC, P c ) and we have shown that this 
statistic estimated on MRLens filtered maps provide a 
strong discrimination between the 5 selected models. 

In this paper, we have compared systematically a wide 
range of non-Gaussian statistics proposed in the literature 
and WPC have been found to be the best statistic. The 
comparison being done on models selected along the (as- 
f2 m )-degeneracy, this study shows that the power spec- 
trum (or the two-point correlation function) and WPC 
should therefore be used simultaneously. However, these 
two statistics probe the same field, therefore the corre- 
lations introduced by the combined measurement n eed to 
be taken into account fsee lTakada and Bridld . I2007T) . This 
will be investigated in the future. 

Another issue, discussed in ^6[ is the selection effects. 
The selection function is needed for the estimation of the 
cosmological parameters and an accurate estimation of 
their errors. This study will be done in a future work. 

Finally, while peak counting and WPC provide a lot of 
information, further statistics such as the cluster count per 
mass, the spatial cluster correlation will provide further 
constraints. Future work will be needed to fully exploit 
this approach. 

Appendix A: The FDR method 

False Discovery Rate (FDR) is a statistical ap- 
pro ach to the multiple testing probl em, introduced 
byfeenj amini and Hochbergl (|1995l ). 

The FDR method offers an effective way to se- 
lect an adaptative threshold, without any assump- 
tion. The FDR threshold is determined from the 
observed p-value distribution, and hence is adap- 
tive to the amount of signal in the data. 

This t echnique has been descr ib ed b y 



Miller et al. 



Starck et al. 



(l200lh : iHopkins et al.1 
(2006]); IPires et ah 



(|2002h : 
(|2006h with 



several examples of astrophysical applications. 
Instead of controlling the chance of any false 
positives, FDR controls the expected proportion 
of false positives. The FDR is given by the ratio 
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(|12p. that is, the proportion of declared active 
which are false positives: 



D a 



(12) 



where Vi a is the number of pixels truly inactive 
declared active, and D a is the number of pixels 
declared active. The FDR formalism ensures that, 
on average, the False Discovery Rate is no larger 
than a which lies between and 1. This procedure 
guarantees control over the FDR in the sense that: 



£{TTJR) < y.a < a 



(13) 



The unknown factor -p- is the proportion of 
truly inactive pixels where Ti is the number of in- 
active pixels and V the total number of pixels. 

The FDR procedure is as follows : 
Let Pi, ...,-P„ denote the p_values from the N tests, 
listed from smallest to largest. 
Let : 



7 Cl 7-> 1 

a = max\ k : Pk < 77 t 

c N .N J 



(14) 



where cn = 1, if p_values are statistically indepen- 
dants. 

Now, declare actived all the pixels with p_values 
less than or equal to Pd. 

Graphically, this procedure corresponds to plot- 
ting the Pfc versus ■£>, superposing the line through 
the origin of slope (see Fig. [9]), and finding the 
last point at which Pf. falls below the line, termed 
Pd- From this p_value Pd, we can derive a threshold 
T. All the pixels greater than T have a p_value less 
than Pd and are declared actives. 



Appendix B: The Higher Criticism definition 

To define HC, first we convert the individual m into p- 
values. Let pi = P{N(0, 1) > Ki} be the ith p-value, and 
let pu) denoted the p-values sorted in increasing order. 
The higher criticism statistic is defined as: 



hc: = 



\fn[i/n-p (l) ] 



y/p^il - P {1) ) 
Or in a modified form : 



(15) 




Figure 9. Finding a threshold graphically using the FDR 
procedure 



HC^ = max 

2:1/ n^P(i) ^ 1 — 1/n 



Vn[i/ n -P(i)} 



(16) 



Appendix C: Description of the representations 

The anisotropic bi-orthogonal wavelet transform 

The most commonly used wavelet transform is the un- 
decimated bi-orthogonal wavelet transform (OWT). Using 
the OWT, a 2D signal S can be decomposed as follows : 



S(x,y) = ^2 <t>J,k x ,k v {x,y)Cj(kx,k y ) 



h h 



EEE 

d k x ,k y j—1 



(17) 



(18) 



with where <j> and ip d are respectively the scaling function 
and the wavelet functions that prioritize the horizontal, 
vertical and diagonal directions. J is the number of res- 
olutions used in the decomposition, Wj the wavelet (or 
detail) coefficients at scale j and direction d, and Cj is a 
smooth version of the original signal S. 

The undecimated isotropic wavelet transform 

The Undecimated Isotropic Wavelet Transform 
(UIWT) decomposes annXtj image k as in : 

,/ 

K(x,y) = Cj(x,y) +22wj(x,y), 
3=1 
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where Cj is a coarse or smooth version of the original 
image k and Wj represents the details in n at scale 2 3 (see 
Starck and Murtagh , 20061 for details). 



The ridgelet transform 

The classical multiresolution ideas only address a por- 
tion of the whole range of interesting phenomena: the 
roughly isotropic one at all scales and all locations. The 
ridgelet transform have been proposed as an alternative 
to the wavelet representation of image data. 

Given a function f{x\,X2), the ridgelet trans- 
form is the superposition of elements of the form 
a _1 / 2 V'((a;i cos 9 + sin# — b)/a), ip is the wavelet, a > 
the scale parameter, b the location parameter and 9 the 
orientation parameter. The ridgelet is constant along lines 
x\ cos 9 + X2 sin 9 — const, and transverse to these ridges 
it is a wavelet. 



The curvelet transform 

Ridgelets are essentially focused on dealing with 
straight lines rather than curves, ridgelets can be adapted 
to representing objects with curved edges using an ap- 
propriate multiscale localization. If one uses a sufficiently 
fine scale to capture curved edges, such edges are almost 
straight. As a consequence the curvelet transform has been 
introduced, in which ridgelet are used in a localized man- 
ner. 



The idea o f th e curvele t 

( Candes and Donohol . 19991 : Starck et al 



tra nsform 
2003h 



is to 

first decompose the image into a set of wavelet planes, 
then to decompose each plane in several blocks (the 
block size can change at each scale level) and to analyse 
each block with a ridgelet transform. The finest the 
scale is, the more sensitive to the curvature the analysis 
is. As a consequence, curved singularities can be well 
approximated with very few coefficients. 



Appendix D: The MRLens filtering 



The MRLens filtering (jStarck et al.l . l2006h 



non linear filtering based on the Bayesian theory 
that searches for a solution that maximizes the a 
posteriori probability. 



Choosing the prior is one of the most critical 
aspects of the Bayesian analysis. The MRLens fil- 
tering uses a multiscale entropy prior. 

Assuming Gaussian noise, the MRlens filtering 
solves the following minimization : 



J(k) 



2ol 



(19) 



3=1 k,l 



where o n the noise standard deviation, J the num- 
ber of scales, f3 is the regularization parameter and 
W is the Wavelet Transform operator. 

Full det ails of the minimiza tion algorithm can 
be found in St arck et al.l ([2001), as well as the way 
to determine automatically the regularization pa- 
rame ter j3. 

In Starck et al. (2006), it has been shown that 
the MRLens filtering outperforms the existing 
methods. The MRLens filtering has already been 
used for several applications on weak lensing 
data and especially, it has been selected to fil- 
ter the dark matter mass map obtained by the 
Hubble Space Telescope in the COSMOS field. 
The complete MRLens software package to per- 
form weak lensing filtering can be downloaded 
from http:/ / www-irfu.cea.fr/Ast/878.html, 
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